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Abstract. The Birnbaum-Saunders distribution, also known as the fatigue-life 
distribution, is frequently used in reliability studies. We obtain adjustments to 
the Birnbaum-Saunders profile likelihood function. The modified versions of 
the likelihood function were obtained for both the shape and scale parameters, 
i.e., we take the shape parameter to be of interest and the scale parameter to be of 
nuisance, and then consider the situation in which the interest lies in performing 
inference on the scale parameter with the shape parameter entering the model- 
ing in nuisance fashion. Modified profile maximum likelihood estimators are 
obtained by maximizing the corresponding adjusted likelihood functions. We 
present numerical evidence on the finite sample behavior of the different estima- 
tors and associated likelihood ratio tests. The results favor the adjusted estima- 
tors and tests we propose. A novel aspect of the profile likelihood adjustments 
obtained in this paper is that they yield improved point estimators and tests. 
The two profile likelihood adjustments work well when inference is made on the 
shape parameter, and one of them displays superior behavior when it comes to 
performing hypothesis testing inference on the scale parameter. Two empirical 
applications are briefly presented. 



1. Introduction 

Birnbaum and Saunders (1969a) derived a two-parameter distribution using a 
set-up in which failure time due to fatigue under cyclic loading when failure fol- 
lows from the development and growth of a dominant crack. According to Mar- 
shall and Olkin (2007), the Birnbaum-Saunders distribution has appeared in sev- 
eral different contexts, and various derivations of the distribution are known. Ac- 
cording to them (pp. 466-467), "it was given by Fletcher (1911), and according 
to Schrodinger (1915) it was obtained by Konstantinowsky (1914);" additionally, 
"it was obtained by Freudental and Shinozuka (1961), but it was the derivation 
of Birnbaum and Saunders (1969a) that brought the usefulness of the distribution 
into clear focus." Desmond (1985) derived the same distribution in a more general 
setting; he used a biological model and relaxed several of the assumptions made 
by the original authors. The relationship between the the Birnbaum-Saunders and 
inverse Gaussian distributions was explored by Desmond (1986). It can shown that 
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the Birnbaum-Saunders distribution is a mixture between an inverse Gaussian dis- 
tribution and a generalized inverse Gaussian distribution; see Bhattacharyya and 
Fries (1982). 

The random variable T is said to be Birnbaum-Saunders distributed, denoted 
T ~ SS(a,/3), if its density function is given by 



t,a,p > 0, where a is the shape parameter and /? is the scale parameter. It is 
noteworthy that the reciprocal property holds for the Birnbaum-Saunders distri- 
bution: ~ SSia,^^); see Saunders (1974). It is easy to show that E(r) = 
/3(l + \a^), var(r) = ia^f(l + fa^), EiT'^) = (l + ^a^) and var(r-i) = 



where <!>(•) denotes the standard normal distribution function. Note that y8 is the 
median of the distribution: FjiP) = 0(0) = 0.5. It was shown by Kundu, Kannan 
and Balakishnan (2008) that the Birnbaum-Saunders hazard function is an upside 
down function for all values of the shape and scale parameters. Hence, the dis- 
tribution is useful in a number of practical situations where the hazard function 
increases up to a point and then decreases. The authors have also addressed the 
important issue of performing inference on the point at which the hazard function 
reaches its maximum. 

Oftentimes the interest lies in performing inference on a subset of the parameters 
that index the model; such parameters are said to be of interest, and the remain- 
ing ones are nuisance parameters. For instance, in Birnbaum-Saunders reliability 
studies, one is typically interested in performing inference on one of the parameters 
that index the model, the other parameter entering the modeling process in nuisance 
fashion. In the presence of nuisance parameters, inferences are usually based on 
the profile likelihood function, which is obtained by replacing, in the likelihood 
function, the nuisances parameters by their corresponding maximum likelihood 
estimators for fixed values of the parameters of interest. The resulting function — 
the profile likelihood function — will only depend on the parameters of interest. 
It is noteworthy, however, that such a function is not a true likelihood function, 
and some of the properties that hold for likelihood functions are no longer valid; in 
particular, there exist score and information biases that do not vanish as the sam- 
ple size increases. Several adjustments to the profile likelihood function have been 
proposed; see, e.g., Barndorff-Nielsen (1983, 1994), Cox and Reid (1987, 1992), 
McCullagh and Tibshirani (1990) and Stern (1997). The main idea behind these 
adjustments is to add a term to the log-likelihood function prior to maximizing it, 
in order to overcome the aforementioned shortcomings. 





The Birnbaum-Saunders distribution function is 
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In this paper we shall use the results in BarndorfF-Nielsen (1983), Severini 
(1998, 1999) and Cox and Reid (1987) to obtain adjustments to the Birnbaum- 
Saunders profile likelihood function. A novel aspect of this approach is that the 
effect of the nuisance parameter on the inference performed on the interest pa- 
rameter is greatly reduced. It is also noteworthy that likelihood ratio type tests 
constructed using the adjusted profile likelihood function typically have superior 
finite sample performance. In short, by adjusting the profile likelihood function 
and then maximizing it one can perform reliable point estimation and hypothesis 
testing inference even when the sample size is small. Our results shall allow prac- 
titioners to perform reliable inference when using the Birnbaum-Sanders model 
in small samples. A motivation for our analysis lies in the important situation in 
which one wishes to make inferences on the the median failure time in a reliability 
study. As we have seen, the median of the Bimbaum-Saunders distribution is /?, 
one of the parameters that index such a distribution. Therefore, here a is a nuisance 
parameter. There are also situations where the interest lies in performing statisti- 
cal inference on the shape parameter a, with yS figuring as a nuisance unknown 
quantity. It is thus important to develop reliable and accurate inference strategies 
that are not sensitive (or, at least, less sensitive) to the parameter that enters the 
modeling in nuisance fashion. This is our chief goal. 

The paper unfolds as follows. Section |2] introduces adjustments to the profile 
likelihood function when the interest lies in performing inference in the presence 
of nuisance parameters. In Section [3l we derive adjustments to the Bimbaum- 
Saunders profile likelihood function. The use of such adjustments delivers, as noted 
above, improved estimation and testing inference in small samples. Alternative 
inference strategies are presented in Section H] Numerical results are presented in 
Sections |5] and m and two applications are presented in Section|7] Finally, Section 
|8] summarizes the main findings and lists directions for future research. 

2. Profile likelihood function and adjustments 

Let , . . . , f„ be independent and identically distributed random variables with 
joint density f{t; 0), where c RP is a /7-vector of unknown parameters and t - 
{t\, . . . , tnY . In what follows, we shall partition G a.?, 9 - (t^, (p^Y , where t, a 
(^-vector, contains the parameters of interest and 4>, 3. ip - ^)-vector, contains the 
nuisance parameters. 

Inference can be based on the profile likelihood function, defined as Lp(T) = 
L(t, (f>r), where L(-) is the usual likelihood function and 4>t is the restricted maxi- 
mum likelihood estimator of 4> given r. The profile likelihood is not a true likeli- 
hood, and some of the properties that hold for a genuine likelihood do not hold for 
its profiled version. In particular, there exist score and information biases, both of 
order 0(1). 

The interest lies in testing the null hypothesis 'Hq : t = tq against 'Hi : t + tq, 
where tq is a given ^-vector of scalars. The likelihood ratio statistic obtained from 
the profile likelihood function is 

LR = 2 {^(?, 0) - €{T,'^r)] = 2 [ip(j) - {p{T)\ . 
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Here, T and ^ are the maximum likelihood estimators of r and (p, respectively, €{■) 
is the log-likelihood function and (p{-) is the profile log-UkeUhood function. Under 
the null hypothesis, LR x^cp where denotes convergence in distribution. 

Several adjustments to the profile likelihood function have been proposed in the 
literature; see, e.g., Severini (2000, Chapter 9), Pace and Salvan (1997, Chapter 
11) and the referenced therein for details. 

Barndorff-Nielsen (1983) proposed an adjusted profile likelihood function which 
is invariant under reparameterizations of the form (t, (p) (t, ^(t, 0)), where t is 
the vector of parameters of interest, (p is the vector of nuisance parameters and ^ is 
a function of r and cp. His proposal follows from the p* formula, which is an ap- 
proximation to the conditional density of the maximum likelihood estimator given 
an ancillary statistic. The proposed adjusted profile likelihood function is 



LBN(r) 



dcPr'' 



d(p 



r,(Pr)\-"%(T), 



where dtprldcp is the matrix of partial derivatives of (pT with respect to (p, j^^(t, (p) - 
-d^€ld(pd(p^ is the observed information matrix for cp when r is fixed and Lp{T) is 
the profile likelihood function for r. 

There is an alternative expression for Lbn that does not involve \d(prld(p\\ it 
involves, nonetheless, a sample space derivative and requires an ancillary statistic 
a such that (r, (p, a) is minimal sufficient. It can be shown that 

— — _j — ^ — 

— ~ = (pT, T, (p, a) i^r^ir, (pr, r, <p, a). 



where 



i^,-^r,<Pr;r,4>,a) = -^\^ J. 



d(p^ 

Here, ^^.^(r, (Pt',t^, (p, a) and j^,p{T, (pj;'T, (p, a) are the log-likelihood function and the 
observed information for <p, respectively. They depend on the data only through the 
minimal sufficient statistic. 

Some approximations to the sample space derivative of the log-likelihood func- 
tion have been proposed. Severini (1998) obtained an approximation to Barndorff- 
Nielsen's adjusted profile UkeUhood function that requires neither a sample space 
derivative nor an ancillary statistic. It is given by 

?Bn{t) ^ {p{T) + ^ log |;^^(t, <pr)\ - log |/^(t, (pr,'?, (p)\, 

where 

/^(t, (p; To, ^o) = E(ro,0o) {U^'^' 4>y^(Jo, (poY] (1) 

is the covariance matrix of log-likelihood derivatives and ^^(t, (p) = d£/d(p. The 
approximation error is of order 0{n~^^^). The corresponding maximum likelihood 
estimator shall be denoted as tbn- 
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An alternative approximation, with the same approximation error, was proposed 
by Severini (1999): 



1 _ „ _ 

^Bw(t) ^ ip{T) + - log |;^0(t, 4>r)\ - log |7^(t, (pr; T, 0)1, 



where 



/^(t, (/>; TO, 0o) = ^J^^^' '^^^"'^(^O' -^o)^' 

7=1 

£fie) = {{[j\e)j'-J\e) being the score function for the jth observation. This ap- 
proximation can be easily computed and is particularly useful in situations where 
one is not able to compute expected values of log-Ukelihood derivatives. The cor- 
responding maximum likelihood estimator shall be denoted as tbn- 

Cox and Reid (1987) defined an adjusted profile likelihood function, where an 
adjustment term is included into the UkeUhood function prior to maximization. It 
approximates the conditional density function of the observations given the nui- 
sance parameter maximum likelihood estimator and can be written as 

Lcr(t) = \U^iTM-"^Lp(T). 

Taking logs we obtain 

^cr(t) = ^(t, 0^) - ^ log \U^(tM (3) 

Note that this function is the penalized counterpart of the log-likelihood function, 

the penalty term taking into account information on the nuisance parameter. The 
maximizer of (cr{t) shall be denoted as tcr. It is noteworthy that the score bias is 
of order 0{rr^), but the information bias remains 0{\). 

The derivation of {cr(t) requires that t and be orthogonal, i.e., that the el- 
ements of the score vector, di/dr and d£/d(f), be uncorrelated which implies that 
iTtf, - 0. When ij^ t 0, it is necessary to find a reparameterization of the form 
(r, A(t, <f>)), where r and A are orthogonal. It is noteworthy that such a reparame- 
terization cannot always be found, except when the parameter of interest is scalar. 
We also note that the Cox and Reid adjustment is not invariant under reparameter- 
izations of the form (r, (f>) (r, ^(r, 0)), unUke Bamdorff-Nielsen's adjustment. 

3. The Birnbaum-Saunders adjusted profile likelihoods 

At the outset, let a be the parameter of interest and /3 the nuisance parameter. 
Also, let t = (fi, . . . , tnV denote a random sample of size n from the Birnbaum- 
Saunders distribution. The log-likelihood function is 

,1/2 /^x3/2l . n , 

_ 1 y + ^ 



£{a,/3) = -n log(a;8) + |] log + 



2 



For fixed a, the restricted maximum likelihood estimator of jS, jSq., is the positive 
root of the following nonUnear equation: 
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where 



1 " 

= - y ti, 



i=l 




and Kifi) = 



i=\ 



-1 



Note that Pa does not have a closed-form expression, and, as a result, it must be 
obtained using restricted nonlinear optimization methods; see, e.g., Nocedal and 
Wright (1999, Chapter 18). (We note that the maximum likelihood estimator of 
yS for fixed a equals the maximum likelihood estimator of that is, ySo, - yS.) By 
replacing (3 hy Pa in i{a,P) we obtain the profile log-likelihood function given by 



i=\ 



J3a 
ti ) 



1/2 



Pa 



3/2 



n 
2a^ 



s 



J3a 



The asymptotic distribution of the vector of maximum likelihood estimators of 
the parameters that index the Birnbaum-Saunders distribution was obtained by 
Englehardt, Bain and Wright (1981). A simpler expression for Fisher's information 
matrix was obtained by Lemonte, Cribari-Neto and Vasconcellos (2007). 

In what follows, we shall obtain the adjusted profile likelihoods described in 
Section |2l Note that the interest and nuisance parameters are orthogonal. The 
adjusted profile log-likelihood function of Cox and Reid (1987) for a can be ex- 
pressed as 

1 



£cR(a) = ^p(a) - - log \jf}/3(a,/3a)\, 



where 



and 



j/s/sia,Pa) 



n 

K'ifi)-- 



+ 



n r 



'acR is the adjusted profile maximum likelihood estimator of a; it does not have 
closed-form and must be obtained numerically. 

The Barndorff-Nielsen (1983) adjusted profile log-likelihood function for a can 
be written as 



hNia) = €p{a) + log 



\j^^{a,l3a)\'l^ 
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Instead of obtaining the term jp.^a, j3a) in £bn^ we shall obtain I{a,]3^; a,P) given 
in (2) using, thus, Severini's (1999) approximation. After some algebra, we obtain 



P P j=i 



1 1 \ f " 2 " 



where 



-3/2 



and 



The adjusted profile maximum likelihood estimator aBN of a cannot be expressed 
in closed-form; it has to computed by numerically maximizing the associated log- 
Ukelihood function. 

The likelihood ratio test statistics obtained from the adjusted profile log-likelihood 
functions given above for the test of "Ho : a - ao against 'Hi : a ^ oq aie 



and 



LRcRia) = 2 {^cRiacR) - ^CR(ao)} 



LRsNia) = 2 [iBN(aBN) - ^Biv(«o)) , 



where Txcr and TiBN are the values of a that maximize (cr{(^) and {bn{o:), respec- 
tively. These test statistics are asymptotically distributed as under the null hy- 
pothesis. 

We shall now consider jS as the parameter of interest and view a as a nuisance 
parameter. For fixed ^, we write the restricted maximum likelihood estimator of a 
as 

By plugging a/j into the log-UkeUhood function we obtain the following profile 
log-likelihood function: 

£p(P) = ^(a^,;S)--^log|^ + ^-2j-nlog;S 

n 

+ J] log 



1=1 



1/2 



.1^ 

ti 



3/2 



The jaa{a,P) block of the observed information matrix evaluated at can 
be written as 

iaa(Si3,P) = -2n 1^+^-2 
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From (3) it follows that Cox and Reid's adjusted profile log-likelihood function for 
jSis 



P s 



The estimator jScR, like the previous estimators, does not have closed-form. 

The Bamdorff-Nielsen adjusted profile log-likelihood function can be expressed 

as 

We use Severini's (1998) approximation and replace Ca\s(!Sp,P), in CsNiP)^ by 
I(afi,fi;a,^) given in (1). We arrive at 



The corresponding estimator, y6g^, does not have closed-form. 

The likelihood ratio test statistics obtained from the above adjusted profile log- 
UkeUhood functions for the test of 'Hq : p - Pq against 9{i : ^ ^ fio axe 



and 



where /3cr and /Sbn are the values of /? that maximize (cRiJ^) and isNiJ^)^ respec- 
tively. The two test statistics are asymptotically distributed d&x\ under "T/q- 



4. Alternative inference strategies 

Some alternative point estimators for the parameters that index the Bimbaum- 
Saunders distributions have been proposed in the literature. Ng, Kundu and Bal- 
akrishnan (2003) obtained modified moment estimators for a and jS. As before, 

let J = i = ^1=1 U (sample arithmetic mean) and r = {n^^ YIi=i (sample 
harmonic mean). The estimators can then be written as 




Ng, Kundu and Balakrishnan (2003) have also proposed jackknife estimators for 
a and yS. The underlying idea is to remove observation tj from the random sample 
t = (fi, f2, . . . , tnY , and to estimate the parameters based on the remaining n - 1 
observations; this is done for j = We shall denote the jackknife maximum 

likelihood estimators as a^gjMLE and ^nsjmle > ^^e jackkinife moment estimator are 

^NgjMME ^^^NgjMME- 



INFERENCE FOR THE BIRNBAUM-SAUNDERS DISTRIBUTION 



9 



From and Li (2006) also proposed alternative estimators for the two parameters 
that index the Bimbaum-Saunders distribution. For instance, they proposed using 



P,, = ^^^ and a,,= .\^M-2 



r 



The authors have also proposed a second estimator for {a,P). Their proposal is 
to estimate /3 using = median(?i, ...,tn) since /3 equals the median of the 
Bimbaum-Saunders distribution. The corresponding estimator for or is 



-2 + 2 Vl + 5v 



where v - ct'^/^fi, being the sample variance, i.e., Tr^ = (n - 1)"^ Tj'i=i(U - 0^- 
Let . . . , ?(„) denote the order statistics of the sample /i, . . . , /„. The estimator 
for jS is, as above, the sample median. For each solve 

i 

F{t(i); a,PF2) 7 , i-\,...,n. 

n + I 

Let i- 1, . . . ,n, denote the solutions, where 

a{i) = 



with h{t) = t^l^ - r^^^. The estimator is aps = median(a(l), . . . , c?(n)). 

From and Li (2006) have also proposed yet another estimator for (a,^). Let 
< A < 0.5, and let «i - nA + I and = n(l - A), to the nearest integer. The 
proposed estimators are 



The authors suggest using A - 0.05, so that only the middle 90% of order statistics 
are used. 

An alternative hypothesis test was proposed by Lemonte, Cribari-Neto and Vas- 
concellos (2007). They derived a Bartlett-correction factor to the likelihood ratio 
statistic and obtained an approximate test whose error vanishes at a faster rate as 
the sample size increases. Let LR* denote their test statistic. It follows that whereas 
Pr(ZJ? < x) = Pr(Yi < x)+ 0(n-^), the correction yields Pr(LR* < x) = Pr(xl < 
x) + 0{n~^), a clear improvement. [See Cribari-Neto and Cordeiro (1996) for a 
detailed review of Bartlett corrections.] Consider the following null hypotheses: 
(i) -Ho : a = 0.1, (ii) 'Ho : a ^ 0.25, (iii) "Hq: a ^ 0.5, (iv) 'Hq: a = 0.75, (v) 
T/o '• a - 1.0, (v) Ho '■ a - 2.0. The corresponding Bartlett-corrected test statistics 
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ar( 



LR* = 



LR 



, LR* = 



LR 



, LR* = 



LR 



LR* = 



1 +4.3918/« 
LR 



, LR* = 



1 + 3.2537/?! 
LR 



and LR* = 



1 + 3.0414/?! ' 
LR 



1 + 2.5924/?! ' 



1 + 2.0307/?! 



1 -0.0445/?!' 



5. Numerical EVIDENCE 



We shall now present Monte Carlo simulation results on the finite sample behav- 
ior of inference based on profile and adjusted profile likelihoods. All simulation 
experiments entail 10,000 replications. The shape parameter assumed two values, 
namely a = 0.5, 1.0, and the scale parameter was set at p = 1.0. The simula- 
tions were performed using the Ox matrix programming language (Doornik, 2006). 
Likelihood maximizations were performed using the quasi-Newton method known 
as BFGS with analytical first derivatives; see Nocedal and Wright (1999) for details 
on the BFGS method. 

Point estimation is evaluated using the following measures: mean, bias, vari- 
ance, mean squared error (MSE), relative bias (RB), asymmetiy and kurtosis. Rel- 
ative bias is defined as 100 x (bias /true parameter value). Hypothesis testing in- 
ference on the parameter of interest is described through the null rejection rates of 
the profile and adjusted profile likelihood ratio tests. Power simulations were also 
performed. 

Table 1 contains simulation results for the case where a is the parameter of 
interest. The sample size is ?! = 10. Note that the estimators c?cR and "Sbn are 
less biased than "Sp. For instance, when a = 0.5 the relative biases ofSp, "Scr and 
c?B/v are 7.50%, 2.16% and 2.13%, respectively. Nevertheless, bias reduction is 
achieved at the expense of greater variability. It is also noteworthy that the small 
sample behavior of the two adjusted profile maximum likelihood estimators are 
similar. We also note that the skewness and kurtosis of Op are slightly closer to 
their asymptotic counterparts than those of 'Scr and TSbn- (When the parameter 
of interest is yS, the estimators /icR and /3bn coincide, since maximization of 
of the profile likelihood function is equivalent to that of icRiP) or {sNifi)- As a 
consequence, the profile and adjusted profile maximum Ukelihood estimators also 
coincide.) 

Table 2 presents the null rejection rates (%) of the different likelihood ratio tests, 
i.e., the tests based on the statistics LR, LRcr, LRbn and LR*, for the test of 'Hq : 
a = ao against "Ki : a ao, where ao is a given scalar; here, a is the parameter 
of interest and p is the nuisance parameter, whose value is set at p = 1.0. The 
values set at the null hypothesis are ao = 0.5, 1.0,2.0 and the sample sizes are 
?! = 10,25,50. All entries are percentages. The figures in Table 2 show that 
the adjusted profile likelihood ratio tests (LRcr and LRbn) outperform the usual 
profile likelihood ratio test (LR). For instance, at the 10% nominal level and when 
oq = 0.5, the rejection rate of the latter was 13.57% whereas the rejection rates 

^Note that the values of the Bartlett correction factor we give correct those given by the authors 
for cases (i) and (ii). 



Table 1 . Point estimation of a. 









a — 


0.5 








estimator 


mean 


bias 


variance 


MSE 


RB(%) 


asymmetry 


kurtosis 


CUp 


0.4625 


-0.0375 


0.0121 


0.0135 


7.5074 


0.2194 


6.0405 


acR 


0.4892 


-0.0109 


0.0136 


0.0138 


2.1695 


0.2355 


6.3629 


aBN 


0.4893 


-0.0107 


0.0137 


0.0138 


2.1385 


0.2356 


6.3648 








a = 


1.0 








estimator 


mean 


bias 


variance 


MSE 


RB(%) 


asymmetry 


kurtosis 


ap 


0.9160 


-0.0840 


0.0472 


0.0543 


8.4020 


0.3715 


11.9221 


acR 


0.9752 


-0.0248 


0.0548 


0.0554 


2.4841 


0.3732 


12.6356 




0.9792 


-0.0208 


0.0567 


0.0572 


2.0834 


0.3733 


12.6829 
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of LRcR and LRbn were 10.85% and 10.86%, respectively. The likelihood ratio 
test is clearly liberal, the null hypothesis being rejected more often than expected. 
The adjusted tests display much smaller size distortions than the likelihood ratio 
test. It is also noteworthy that the finite sample behavior of the likelihood ratio 
test deteriorates as the value of the shape parameter increases, especially when the 
sample size is small; the adjusted tests remain reliable. The Bartlett-corrected test 
is clearly outperformed by the adjusted profile likelihood tests when a is large and n 
is small. For example, when a = 2.0, n = 10 and the nominal level is 5%, the null 
rejection rate of the Bartlett-corrected test is 9.40% whereas the adjusted profile 
likelihood tests reject the null hypothesis 5.20% (LRcr) and 4.58% (LRbn) of the 
time. Note also that the tests based on LRcr and LRbn display similar small sample 
behavior, especially when the value of the shape parameter is small. (Values of a 
between 0.1 and 0.5 are common in fatigue studies.) 

We have also performed simulations under the alternative hypothesis. The pow- 
ers of the tests LR, LRcr, LRbn and LR* at the 5% and 1% nominal levels were 
computed for values of a ranging from 0.12 to 0.28, the null hypothesis under test 
being 'Hq : a = a^. The tests were carried out using exact critical values, which 
were estimated in the size simulations. This was done so that the different tests 
have the same size, and power comparisons become meaningful. The results are 
presented in Table 3 and were obtained using n = 10, a - 0.10 and ^ - 1.0. (All 
entries are percentages.) We note that LR is slightly less powerful than LRbn and 
LRcr. For example, when a - 0.20 and at the 5% nominal level, the nonnuU rejec- 
tion rates of these tests were equal to 77.86%, 83.81% and 83.81%, respectively; 
the nonnull rejection rate of the Bartlett-corrected test was 77.87%. 

Figure 1 plots the relative quantile discrepancies of the three test statistics against 
the corresponding asymptotic quantiles. Relative quantile discrepancy is defined 
as the difference between exact (estimated by simulation) and asymptotic quan- 
tiles divided by the latter. The closer to zero the relative quantile discrepancies, the 
better the approximation of the exact null distribution of the test statistic by the lim- 
iting distribution. It is noteworthy that the relative quantile discrepancies of the 
adjusted test statistics are considerably closer to zero than those of the likelihood 
ratio test statistic, which oscillate around 18%. The relative quantile discrepancies 
of the two adjusted test statistics are very similar. 

Figure 2 plots the relative size distortions against the corresponding nominal 
levels of the tests. Relative size distortion is defined as the difference between p- 
values (estimated by simulation) and nominal levels divided by the latter. Note 
that the relative size distortion of the likelihood ratio test increases rapidly as the 
nominal level of the test decreases, which does not occur for the adjusted tests. 
Note also that the relative size distortions of the two adjusted tests are very similar. 

Table 4 contains the null rejection rates (again expressed as percentages) of the 
three tests (LR, LRcr and LRbn) for the test 'Ho '■ P - Po- (Note that here p is 
the parameter of interest.) Since fi functions as a multiplier, as explained earlier, 
we have only performed simulations using yS = 1 .0; four different values of a were 
used, namely: 0.1, 0.5, 1.0 and 2.0. As in Table 2, three different sample sizes 
were considered: n - 10,25,50. The figures in Table 4 reveal that the likelihood 



Table 2. Null rejection rates, inference on a (J3 - 1.0). 





nominal 




a = 


0.1 






a = 


0.5 






a = 


1.0 






a = 


2.0 




n 


level 


LR 


LRcR 


LRbn 


LR* 


LR 


LRcR 


LRbn 


LR' 


LR 


LRcR 


LRbn 


LR' 


LR 


LRcR 


LRbn 


LR* 




10 


13.44 


11.00 


11.00 


7.24 


13.57 


10.85 


10.86 


9.03 


13.97 


10.65 


10.78 


10.61 


16.10 


9.89 


8.99 


15.90 


10 


5 


7.47 


5.14 


5.14 


3.05 


7.60 


5.18 


5.20 


3.96 


7.83 


5.24 


5.29 


5.11 


9.44 


5.20 


4.58 


9.40 




1 


1.70 


1.07 


1.07 


0.45 


1.72 


1.11 


1.11 


0.69 


1.83 


1.06 


1.11 


0.93 


2.53 


1.04 


1.96 


2.51 




0.5 


0.89 


0.53 


0.53 


0.16 


0.91 


0.52 


0.52 


0.31 


0.97 


0.52 


0.55 


0.51 


1.40 


0.48 


0.40 


1.39 




10 


10.84 


10.15 


10.15 


8.42 


10.88 


10.13 


10.14 


9.15 


11.56 


10.41 


10.50 


10.30 


11.98 


10.12 


10.10 


11.96 


25 


5 


5.94 


5.33 


5.33 


4.32 


5.93 


5.34 


5.35 


4.67 


6.18 


5.31 


5.33 


5.28 


6.24 


5.36 


5.30 


6.22 




1 


1.43 


1.18 


1.18 


0.83 


1.46 


1.15 


1.15 


1.02 


1.46 


1.04 


1.04 


1.04 


1.54 


1.07 


1.06 


1.54 




0.5 


0.82 


0.57 


0.57 


0.43 


0.82 


0.59 


0.59 


0.48 


0.61 


0.50 


0.52 


0.43 


0.82 


0.59 


0.54 


0.81 




10 


10.89 


10.35 


10.35 


9.46 


10.86 


10.35 


10.34 


9.99 


10.93 


10.31 


10.32 


10.25 


11.24 


10.48 


10.47 


11.23 


50 


5 


5.54 


5.21 


5.21 


4.58 


5.53 


5.20 


5.19 


4.81 


5.40 


5.03 


5.05 


4.86 


5.65 


4.95 


5.04 


5.64 




1 


1.19 


0.99 


0.99 


0.87 


1.21 


1.02 


1.02 


0.99 


1.14 


1.02 


1.02 


1.10 


1.18 


1.03 


1.03 


1.17 




0.5 


0.70 


0.56 


0.56 


0.51 


0.71 


0.55 


0.55 


0.58 


0.68 


0.61 


0.61 


0.61 


0.60 


0.49 


0.51 


0.60 
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nominal level: 5% 



nominal level: 1 % 



a 



BN 



LR* 



LR LR(^]i LR 



BN 



LR* 



0.12 
0.14 
0.16 
0.18 
0.20 
0.22 
0.24 
0.26 
0.28 



8.55 
24.79 
45.74 
64.78 
77.86 
86.41 
91.25 
94.91 
96.98 



13.13 
34.07 
56.06 
72.63 
83.81 
89.85 
94.18 
96.73 
97.95 



13.13 
34.07 
56.06 
72.63 
83.81 
89.85 
94.18 
96.73 
97.95 



8.56 
24.79 
45.75 
64.58 
77.87 
86.41 
91.25 
94.91 



2.50 
11.71 
29.37 
48.57 
65.33 
77.29 
85.55 
90.17 
93.90 



4.87 
18.01 
38.16 
59.91 
72.58 
82.97 
88.88 
93.01 
95.40 



4.87 

18.01 

38.16 

59.91 

72.58 

82.97 

88.88 

93.01 

95.40 



2.50 
11.71 
29.37 
48.57 
65.33 
77.29 
85.55 
90.17 
93.90 



ratio test is liberal, that LRbn is conservative, and that LRcr displays very minor 
size distortions, the latter clearly outperforming the other tests. For instance, when 
n = 10, a = 2.0 and at the 10% nominal level, the null rejection rates of these tests 
are, respectively, 13.12%, 7.94% and 10.66%. 

In Table 5 we present the empirical powers of the tests of 'Ho '■ P - Po- The 
values of /J used ranged from 1.2 to 4.0. Again, the tests were performed using 
size-corrected critical values (obtained from the size simulations) in order to force 
them to have the correct size. The simulations were carried out using « = 10, 
Of = 1 .0 and ^ - 1.0. The results suggest that the powers of the thi^ee tests are very 
similar, with a slight advantage of LR. For example, when /J - 2.0, the nonnull 
rejection rates of LR, LRqr and LRbn at the 5% nominal level are, respectively, 
58.54%, 58.05% and 57.43%. 

Figure 3 plots the relative quantile discrepancies against the corresponding as- 
ymptotic quantiles, and Figure 4 plots the relative size distortions against the cor- 
responding nominal levels of the tests when (5 is the parameter of interest. Both 
figures clearly show that LRcr outperforms LR and LRbn- 

6. Additional numerical evidence: comparison with alternative estimators 

We shall now compare the small sample behavior of our adjusted profile max- 
imum likelihood estimators of a to those proposed by Ng, Kundu and Balakrish- 
nan (2003) and From and Li (2006), which are described in SectionUl The number 
of Monte Carlo replications is, as before, 10,000, the true values of a are 0.5 and 
1.0, and n = 10. The simulation results are presented in Tabled 

The figures in Table [6] show that no estimator uniformly outperforms all others 
in terms of both bias and mean squared error. They also show that the adjusted 
profile maximum likelihood estimators are competitive, since they are amogst the 
best performing estimators in both situations (a = 0.5 and a = 1.0). When a = 
0.5, the least biased estimator is af4 (relative bias: 0.66%) whereas when the true 
parameter value is 1.0 the estimator ap^ has the smallest relative bias (0.26%). 
In both cases, "Sbn is the estimator with the third smallest relative bias, followed 



Table 4. Null rejection rates, inference on j3 (fi - 1.0). 





nominal 


a = 0.1 


or = 0.5 


a= 1.0 


or = 2.0 


n 


level 


LR 


LRcR 


LRbn 


LR 


LRcR 


LRbn 


LR 


LRcR 


LRbn 


LR 


LRcR 


LRbn 




10 


12.31 


10.30 


8.52 


12.33 


10.25 


8.26 


12.37 


10.18 


8.01 


13.12 


10.66 


7.94 


10 


5 


6.61 


5.35 


3.95 


6.62 


5.30 


3.85 


6.75 


5.18 


3.65 


6.95 


5.27 


3.65 




1 


1.49 


1.08 


0.70 


1.56 


1.05 


0.69 


1.55 


1.09 


0.67 


1.73 


1.12 


0.54 




0.5 


0.91 


0.59 


0.35 


0.89 


0.60 


0.33 


0.87 


0.59 


0.28 


0.97 


0.53 


0.25 




10 


11.10 


10.38 


9.67 


11.19 


10.33 


9.60 


11.15 


10.37 


9.52 


11.27 


10.52 


9.53 


25 


5 


5.83 


5.30 


4.86 


5.91 


5.38 


4.76 


5.87 


5.33 


4.75 


6.14 


5.49 


4.77 




1 


1.25 


1.08 


0.93 


1.30 


1.07 


0.99 


1.40 


1.15 


0.97 


1.19 


1.09 


0.96 




0.5 


0.62 


0.52 


0.44 


0.64 


0.52 


0.45 


0.76 


0.61 


0.46 


0.80 


0.64 


0.43 




10 


10.71 


10.43 


9.96 


10.64 


10.35 


9.92 


10.72 


10.24 


9.74 


10.87 


10.38 


9.87 


50 


5 


5.32 


5.05 


4.85 


5.36 


5.17 


4.88 


5.54 


5.17 


4.90 


5.39 


5.07 


4.79 




1 


1.19 


1.08 


1.00 


1.20 


1.10 


0.99 


1.18 


1.09 


0.98 


1.12 


1.06 


0.97 




0.5 


0.56 


0.51 


0.48 


0.58 


0.53 


0.46 


0.61 


0.55 


0.50 


0.65 


0.59 


0.52 
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Table 5. NonnuU rejection rates, inference on j6. 





nominal level 


: 5% 


nominal level 


: 1% 


p 


LR 


LRcR 


LRbn 


LR 


LRcR 


LRbn 


1.2 


9.17 


9.07 


9.06 


2.06 


2.06 


2.07 


1.6 


31.67 


31.32 


31.06 


10.32 


10.18 


9.94 


2.0 


58.54 


58.05 


57.43 


24.56 


24.21 


23.50 


2.4 


80.11 


79.72 


79.13 


45.06 


44.07 


42.93 


2.8 


91.05 


90.93 


90.49 


62.15 


61.08 


59.56 


3.2 


96.80 


96.53 


96.42 


75.86 


74.73 


72.68 


3.6 


99.05 


99.01 


98.85 


85.86 


84.74 


83.10 


4.0 


99.61 


99.55 


99.45 


92.34 


91.29 


89.75 



by c?cR- When a - 0.5, for instance, the relative biases of these estimators are 
nearly four times smaller than those of the jackknife estimators and nearly 3.5 
times smaller than the relative bias of the modified moments estimator. 

We note that the approach proposed in this paper, namely adjusting the profile 
log-likelihood function prior to maximization, has a clear advantage over the al- 
ternative approaches described in Section IH it not only improves the small sample 
performance of point estimators, but also improves the finite sample behavior of as- 
sociated likelihood ratio tests. That is, the correction delivers improved estimation 
and testing inference in small samples. 

7. Applications 

We shall now perform profile and adjusted profile Ukelihood inference using two 
real data sets. In both cases, we shall assume that observations are random draws 
from the Birnbaum-Saunders distribution. 

At the outset, we consider the data provided by Birnbaum-Saunders (1969b) 
on the fatigue life of 6061-T6 aluminum coupons cut parallel to the direction of 
rolling and oscillated at 18 cycles per second (cps). The data set consists of 101 
observations with maximum stress per cycle 31,000 psi. Let a be the parameter 
of interest. The profile and adjusted profile maximum likelihood estimates are 
S" = 0.17038, "ScR = 0.17125 md'oBN = 0.17122. Suppose we are interested in 
testing TYo : a = 0.15 against 'Hi : a 0.15. The test statistics based on ip{a), 
^cr{o:) and fsNio:) are, respectively, 3.5771, 3.8421 and 3.8351, with the following 
corresponding jo-values: 0.05858, 0.04998 and 0.05019. Since the sample size is 
large (101 observations), the values of the three statistics are similar. However, the 
resulting inference is not the same at the 5% nominal level, since the test based on 
{cR{a), unlike the other two tests, yields rejection of the null hypothesis. 

We shall now turn to the case where p is the parameter of interest. In particular, 
we are interested in testing "Ko '■ P = 125. The test statistics are LR = 9.4279, 
LRcR = 9.3338 and LRbn - 9.2397, with the following corresponding /^-values: 
0.00214, 0.00225 and 0.00237. Unlike the previous inference, here the three tests 



Table 6. Point estimation of a revisited. 









or - 


0.5 








estimator 


mean 


bias 


variance 


MSE 


RB(%) 


asymmetry 


kurtosis 




0.4625 


-0.0375 


0.0121 


0.0135 


7.5074 


1.3537 


5.9906 


acR 


0.4892 


-0.0109 


0.0136 


0.0138 


2.1695 


1.4258 


6.3098 


aBN 


0.4893 


-0.0107 


0.0137 


0.0138 


2.1385 


1.4624 


6.3117 


0!MME 


0.4625 


-0.0375 


0.0121 


0.0135 


7.5077 


1.3537 


5.9906 


aNg 


0.4749 


-0.0251 


0.0127 


0.0133 


5.0275 


1.3874 


6.1376 


'^NgjMLE 


0.4573 


-0.0427 


0.0133 


0.0152 


8.5384 


1.3397 


5.9301 


^NgjMME 


0.4579 


-0.0421 


0.0133 


0.0151 


8.4255 


1.3412 


5.9367 


api 


0.4625 


-0.0375 


0.0121 


0.0135 


7.4949 


1.3539 


5.9913 


api 


0.4690 


-0.0310 


0.0166 


0.0176 


6.2008 


1.3715 


6.0678 


aps 


0.5065 


0.0065 


0.0326 


0.0326 


1.3004 


1.4721 


6.5222 


apA 


0.5033 


0.0033 


0.0220 


0.0220 


0.6551 


1.4635 


6.4825 








a = 


1.0 








estimator 


mean 


bias 


variance 


MSE 


RB(%) 


asymmetry 


kurtosis 


ap 


0.9160 


-0.0840 


0.0472 


0.0543 


8.4020 


2.3784 


11.8471 


acR 


0.9752 


-0.0248 


0.0548 


0.0554 


2.4841 


2.4786 


12.5612 




0.9792 


-0.0208 


0.0567 


0.0572 


2.0834 


2.4852 


12.6085 


0!MME 


0.9158 


-0.0842 


0.0471 


0.0542 


8.4173 


2.3782 


11.8453 


aNg 


0.9404 


-0.0596 


0.0497 


0.0533 


5.9615 


2.4207 


12.1452 


^NgjMLE 


0.9050 


-0.0950 


0.0523 


0.0613 


9.4959 


2.3591 


11.7120 


^NgjMME 


0.9061 


-0.0939 


0.0522 


0.0611 


9.3931 


2.3609 


11.7247 


&Pl 


0.9168 


-0.0832 


0.0475 


0.0544 


8.3177 


2.3799 


11.8575 




0.8916 


-0.1084 


0.0763 


0.0881 


10.8379 


2.3350 


11.5450 


apj 


0.9974 


-0.0026 


0.1279 


0.1279 


0.2598 


2.5144 


12.8217 


ap^ 


1.0039 


0.0039 


0.0906 


0.0906 


0.3915 


2.5247 


12.8972 
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yield the same conclusion: the null hypothesis is rejected at the 1 % nominal level. 
The (profile and adjusted profile) maximum likelihood estimate of jS is 131.8188. 

The second application we consider uses data provided by McCool (1974). The 
data describe the lifetime, in hours, of 10 sustainers of a certain type. They were 
used by Cohen, Whitten and Ding (1984) as an illustration of the three-parameter 
WeibuU distribution. The profile and adjusted profile maximum likelihood esti- 
mates of (2 are (? = 0.2825, c?cr - 0.2989 and Iibn - 0.2973. Consider the test 
ofHo : a = 0.21 against "Ki : a 0.21. The test statistics are L/? = 2.1646, 
LRcR - 2.8438 and LRbn - 2.7963, with the following corresponding p-values: 
0.1412, 0.0917 and 0.0945. Therefore, the two adjusted profile likelihood ratio 
tests (LRcR and LRbn) reject the null hypothesis at the 10% nominal level, unlike 
the likelihood ratio test (LR). Thus, the unadjusted and adjusted tests yield different 
conclusions. 

Now let jS be the parameter of interest. Its maximum likelihood estimate is 
is = 212.05. Suppose we wish to test 'Ho : /3 = 180 against : P 180. The test 
statistics are LR = 2.9417, LRcr = 2.6415 and LRbn = 2.3414; the respective p- 
values are 0.0863, 0.1041 and 0.1260. Again, the use of adjustments to the profile 
likelihood function makes a difference: unlike the usual likelihood ratio test, the 
adjusted likelihood tests reject the null hypothesis at the 10% nominal level. 

8. Concluding remarks 

We considered the issue of performing inference on the parameters that index 
the Bimbaum-Saunders distribution. More specifically, we have focused on the 
situation where one wishes to make inference on one of the parameters, the other 
parameter being of nuisance fashion. Using the results in Cox and Reid (1987) and 
in Bamdorff-Nielsen (1983), we derived two adjustments that can applied to the 
profile likelihood function so as to deliver improved inference. Approximations 
due to Severini (1998, 1999) were used in order to obtain one of such adjustments. 
Monte Carlo simulation results have shown that the adjusted estimators and tests 
— i.e., estimators and tests based on the adjusted profile likelihood functions — 
can deliver more accurate inference than that carried out using the usual maximum 
likelihood estimator and the standard likelihood ratio test in small samples. In 
particular, the adjusted estimators displayed smaller biases and the adjusted tests, 
smaller size distortions. For instance, we reported Monte Carlo simulation results 
in which the usual likelihood ratio test displayed null rejection rate of nearly 9.5% 
at the 5% nominal level whereas the sizes of our two adjusted tests where 5.2% and 
4.6%, and in which the relative bias of our two estimators were approximately four 
times smaller than that of the maximum likelihood estimator. The adjusted profile 
likelihood tests have also outperformed the Bartlett-corrected likehhood ratio test. 
We recommend the use of the adjusted profile likelihood inference developed in 
this paper to practitioners who wish to model reliabiUty data using the Birnbaum- 
Saunders model. In particular, we recommend the use of the Cox-Reid adjusted 
profile likelihood function, since it yielded the most reliable (hypothesis testing) 
inference when the parameter was of interest was p and it was competitive with the 
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inference obtained using the BarndorfF-Nielsen modified profile likelihood func- 
tion when a was the parameter of interest. In future research, we shall obtain 
adjustments to Birnbaum-Saunders profile likelihoods under data censoring. 
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Figure 2. Relative size distortion plot, inference on a. 
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Figure 3. Relative quantile discrepancy plot, inference on jS. 




Figure 4. Relative size distortion plot, inference on jS. 



